In [1]:
import os
import re
import math
import tifffile
import cv2
import sys
import numpy as np
import matplotlib.pyplot as plt
from collections import Counter
from scipy.stats import scoreatpercentile
from apeer_ometiff_library import io
from pathlib import Path
sys.path.append("..")
import src.ImageParser as IP
In [2]:
# path_images_raw = './data/IMC_ESD/2.ROIs_raw/'
path_images_raw = '../data/IMC_ESD/raw'
fovs = os.listdir(path_images_raw)
print(fovs)
['ESD4_02', 'ESD4_03', 'ESD1_05', 'ESD2_05', 'ESD8_05', 'ESD7_03', 'ESD7_08', 'ESD1_02', 'ESD7_04', 'ESD6_03', 'ESD4_05', 'ESD7_07', 'ESD2_03', 'ESD3_08', 'ESD3_06', 'ESD3_07', 'ESD8_04', 'ESD7_06', 'ESD1_03', 'ESD5_05', 'ESD8_03', 'ESD5_03', 'ESD6_07', 'ESD4_01', 'ESD2_02', 'ESD6_04', 'ESD6_02', 'ESD5_02', 'ESD3_03', 'ESD2_06', 'ESD2_01', 'ESD4_04', 'ESD8_09', 'ESD1_04', 'ESD7_09', 'ESD3_05', 'ESD1_08', 'ESD3_01', 'ESD8_01', 'ESD7_05', 'ESD8_07', 'ESD6_08', 'ESD8_08', 'ESD5_01', 'ESD5_04', 'ESD3_02', 'ESD3_04', 'ESD6_09', 'ESD5_06', 'ESD6_06', 'ESD7_01', 'ESD6_05', 'ESD8_06', 'ESD2_04', 'ESD1_09', 'ESD7_02', 'ESD1_01', 'ESD8_02', 'ESD6_01', 'ESD1_06', 'ESD1_07']
In [3]:
channels = ['CD20', 'Vimentin', 'PD-L1', 'CD31', 'CD163', 'VISTA', 'Ki-67', 'DNA2', 'IDO',
'FOXP3', 'CD68', 'CD57', 'CD14', 'D2-40', 'CD56', 'CD45RO', 'DNA1', 'CD11c', 'CD7',
'HLA-DR', 'CD204', 'CD8a', 'P16Ink4a', 'CD3', 'Granzyme B', 'Bcatenin',
'Caspase', 'CD4', 'CD103', 'TGFbeta', 'PD-1', 'CD45', 'LAG-3',
'ICOS', 'CD11b', 'Keratin', 'TCRgd', 'CD15', 'TIM-3', 'CD38', 'Tbet', 'CD39']
channel_to_index = {name: index for index, name in enumerate(channels)}
def get_files_from_dir(fov):
fov_folder = os.path.join(path_images_raw, fov)
fov_files = os.listdir(fov_folder)
fov_files = [filename for filename in fov_files if filename.lower().endswith((".tiff", ".tif"))]
fov_files = [os.path.join(fov_folder, file) for file in fov_files]
return fov_files
def get_stack(fov_files):
# get a stack of images
stack = []
for channel in channels:
file = [filename for filename in fov_files if channel in filename][0]
img_apeer, _ = io.read_ometiff(file)
img = img_apeer.squeeze()
stack.append(np.array(img))
stack = np.stack(stack, axis = -1)
return stack
def plot_comparison(image1, image2, channel_to_compare):
plt.figure(figsize=(10, 5), dpi = 500)
plt.subplot(1, 2, 1)
plt.imshow(image1, cmap='gray',vmin = 0, vmax = 1) #
plt.title(channel_to_compare)
plt.subplot(1, 2, 2)
plt.imshow(image2, cmap = 'gray', vmin = 0, vmax = 1) # , vmin = 0, vmax = 1
plt.title(channel_to_compare)
plt.show()
def save_stack(img, output_dir):
if not os.path.exists(output_dir):
os.makedirs(output_dir)
for i in range(img.shape[2]):
slice = img[:,:,i]
channel = channels[i]
filename = os.path.join(output_dir, f"{channel}.tiff") # .png
# Save the image
tifffile.imwrite(filename,np.float32(slice),
photometric="minisblack")
# plt.imshow(np.float32(slice), cmap='gray') # Assuming grayscale image, adjust cmap if needed
# plt.axis('off') # Turn off axes
# # Save the image as a PNG file
# plt.savefig(os.path.join(output_dir, f"{channel}.png"), bbox_inches='tight', pad_inches=0,dpi = 700)
# plt.close()
# functions for neigbooring and counting pixels analysis
def get_counts_neigh(image):
# Threshold the image to get binary values
binary_image = (image > 0).astype(np.uint8)
positive_pixel_count = np.count_nonzero(binary_image)
# print(positive_pixel_count)
# Find the coordinates of positive (non-zero) pixels in the binary image
positive_pixel_coords = np.argwhere(binary_image == 1)
# Initialize lists to store results
positive_counts = []
medians = []
percentile_25 = []
percentile_75 = []
# Iterate through positive pixels
for coord in positive_pixel_coords:
y, x = coord[0], coord[1]
# Extract the 3x3 window around the current positive pixel
# Define the coordinates for the 3x3 window
y_start, y_end = max(y - 1, 0), min(y + 2, binary_image.shape[0])
x_start, x_end = max(x - 1, 0), min(x + 2, binary_image.shape[1])
# Extract the 3x3 window around the current positive pixel
window = binary_image[y_start:y_end, x_start:x_end]
# Count positive pixels in the 3x3 window
positive_count = np.count_nonzero(window)
positive_counts.append(positive_count)
# Calculate the median and percentiles
window = image[y_start:y_end, x_start:x_end]
median = np.median(window)
p_25 = scoreatpercentile(window, 25)
p_75 = scoreatpercentile(window, 75)
# Append the statistics to the respective lists
medians.append(median)
percentile_25.append(p_25)
percentile_75.append(p_75)
return positive_counts, medians, percentile_25, percentile_75
def plot_hist_positive_neig(positive_counts1,positive_counts2):
# Create a figure with two subplots (side by side)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# Plot the first histogram on the left subplot
axes[0].hist(positive_counts1, bins=range(max(positive_counts1) + 2), rwidth=0.8, align='left')
axes[0].set_xlabel('Number of neighboring Positive Pixels of positive pixels')
axes[0].set_ylabel('Frequency')
axes[0].set_title('Histogram 1')
# Plot the second histogram on the right subplot
axes[1].hist(positive_counts2, bins=range(max(positive_counts2) + 2), rwidth=0.8, align='left')
axes[1].set_xlabel('Number of neighboring Positive Pixels of positive pixels')
axes[1].set_ylabel('Frequency')
axes[1].set_title('Histogram 2')
# Adjust spacing between subplots
plt.tight_layout()
# Show the plots
plt.show()
def plot_percentiles(list1,list2):
# Create two separate figures for the first and second graphs
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
# Plot the distribution of medians and percentiles for the first graph
ax1.hist(list1, bins=10, color=['blue', 'green', 'red'], alpha=0.7, label=['Median', 'Percentile 25', 'Percentile 75'])
ax1.set_xlabel('Values')
ax1.set_ylabel('Frequency')
ax1.set_title('Distribution of Medians and Percentiles ()')
ax1.legend()
# Plot the distribution of medians and percentiles for the first graph
ax2.hist(list2, bins=10, color=['blue', 'green', 'red'], alpha=0.7, label=['Median', 'Percentile 25', 'Percentile 75'])
ax2.set_xlabel('Values')
ax2.set_ylabel('Frequency')
ax2.set_title('Distribution of Medians and Percentiles ()')
ax1.legend()
# Adjust spacing between subplots
plt.tight_layout()
# Show the plots
plt.show()
In [80]:
fov = 'ESD1_01/'
fov_files = get_files_from_dir(fov)
stack = get_stack(fov_files)
print(stack.shape)
# apply percentile
img_out = IP.remove_outliers(stack, up_limit=99, down_limit=1)
# normalize
img_norm = IP.normalize_by_channel(img_out)
print(np.max(img_norm))
print(np.min(img_norm))
(1000, 1000, 42) 1.0 0.0
In [8]:
def filter_hot_pixelsBodenmiller(img: np.ndarray, thres: float) -> np.ndarray:
'''
:param img:
:param thres:
:return:
# https://github.com/BodenmillerGroup/ImcSegmentationPipeline/blob/56ce18cfa570770eba169c7a3fb02ac492cc6d4b/src/imcsegpipe/utils.py#L10
'''
from scipy.ndimage import maximum_filter
kernel = np.ones((1, 3, 3), dtype=bool)
kernel[0, 1, 1] = False
# array([[[ True, True, True],
# [ True, False, True],
# [ True, True, True]]])
max_neighbor_img = maximum_filter(img, footprint=kernel, mode="mirror")
return np.where(img - max_neighbor_img > thres, max_neighbor_img, img)
In [69]:
img_fil = filter_hot_pixelsBodenmiller(img_norm, thres = 0.05)
channel_to_compare = 'Bcatenin'
index = channel_to_index.get(channel_to_compare, -1)
plot_comparison(stack[:,:,index], img_fil[:,:,index], channel_to_compare)
channel_to_compare = 'CD20'
index = channel_to_index.get(channel_to_compare, -1)
plot_comparison(stack[:,:,index], img_fil[:,:,index], channel_to_compare)
